R Notebook

Step 1: revised catch for 2010-2020 Step 2: tune resilience of all the species (i.e., long-lived, slower growing species = vulnerable…etc…) Step 3: calibrated model with fishing mortality represented (time-averaged catch comparison, model vs observed) = base model Step 4: check base model without fishing to ensure co-existence is still present, observed biomass would not be specified Step 5: Experiments through time - effort functions - uses the base model with no fishing effort as a starting point Step 6: Expose the base model to fishing through time and compare to observed catch time series - start from 0 effort at first time-step to max effort approx around collapse - If the modelled versus observed don’t match, need to adjust effort/or initial biological params and experiment to see what improves - hand tune / optim options Adjust and recalibrate if needed (optim) Use rules to constrain - i.e., - time-averaged biomasses within calibration period must be +- xx % - can optim things like fishing-size selectivity, PPMRs etc

Can we reproduce the time series?

Load libraries

remotes::install_github("sizespectrum/mizerExperimental")

Skipping install of 'mizerExperimental' from a github remote, the SHA1 (8279ac0d) has not changed since last install.
  Use `force = TRUE` to force installation

library(mizerExperimental)
# remotes::install_github("sizespectrum/mizerMR")
# library(mizerMR)
# library(mizer)
library(tidyverse)

Warning: package ‘tidyverse’ was built under R version 4.2.3Warning: package ‘ggplot2’ was built under R version 4.2.3Warning: package ‘readr’ was built under R version 4.2.3Warning: package ‘forcats’ was built under R version 4.2.3── Attaching core tidyverse packages ──────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse 2.0.0 ──
✔ dplyr     1.1.0     ✔ readr     2.1.4
✔ forcats   1.0.0     ✔ stringr   1.5.0
✔ ggplot2   3.4.2     ✔ tibble    3.1.7
✔ lubridate 1.9.2     ✔ tidyr     1.3.0
✔ purrr     1.0.1     ── Conflicts ────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
✖ dplyr::filter() masks stats::filter()
✖ dplyr::lag()    masks stats::lag()
ℹ Use the ]8;;http://conflicted.r-lib.org/conflicted package]8;; to force all conflicts to become errors

library(tictoc)
library(parallel)
# library(plotly)

Catch data for all indivudals and summarised as totals per species each year

Maybe this isn’t expressing effort, rather fishing mortality

# df_ind_CPUE <- readRDS("ind_catch_weight_BanzareBank_1930_2019_CPUE.rds")
# df_CPUE_kg_day <- readRDS("catch_timeseries_BanzareBank_1930_2019_CPUE.rds")
# 
# glimpse(df_ind_CPUE)
# glimpse(df_CPUE_kg_day)

Yield for the period that matches Ecopath model, post 2000 annual average

# df_yield_2000 <- df_CPUE_kg_day %>% 
#   filter(Year > 2009) %>% 
#   group_by(Species) %>% 
#   summarise(yield = sum(total_catch_kg)/10)
# 
# df_yield_2000

effort_days/max value plotted by species This might indicate the effort curve required, especially if there is a difference between different species

# df_plot <- df_CPUE_kg_day %>% 
#   filter(Year > 2000) %>% 
#   group_by(Species) %>% 
#   mutate(max_effort = max(effort_days)) %>% 
#   group_by(Species,Year) %>% 
#   mutate(effort_standard = effort_days/max_effort)

Relative change in catch pre-1930, could inform the relative change in effort Ask Cami, how could I reconstruct the effort to the initial point of catch when effort = 0…i.e., approx 1850

# df_plot %>% 
#   ggplot(aes(x = Year, y = effort_standard, colour = Species)) +
#   geom_smooth() +
#   geom_line()
#   # facet_wrap(~Species) 

# df_plot %>% 
#   ggplot(aes(x = Year, y = CPUE, colour = Species)) +
#   # geom_smooth() +
#   geom_line()
#   # facet_wrap(~Species) 

Load catch length

catch_lengths <- readRDS("catch_lengths.rds")

glimpse(catch_lengths)

The value for euphausiids obs_yield is from calibration_catch_histsoc_1850_2004_regional_models.csv found at http://portal.sf.utas.edu.au/thredds/catalog/gem/fishmip/ISIMIP3a/InputData/fishing/histsoc/catalog.html. It is the annual average yield over a 22 year period for the Prydz Bay Region.

The values of obs_yield for minke whales (2175476136), orca (19923729), sperm whales (4062177445), and baleen whales (50341696055) are the annual average yield from 1930 - 2019 from IWC records of whaling in the Prydz Bay model domain (1,433,028 km2).

Adjust the catch to represent the 2010 - 2020 time period + incorporate toothfish catch if possible (i.e., Stacey’s paper values)

df_ind_CPUE <- readRDS("ind_catch_weight_BanzareBank_1930_2019_CPUE.rds")
df_CPUE_kg_day <- readRDS("catch_timeseries_BanzareBank_1930_2019_CPUE.rds")

glimpse(df_ind_CPUE)
glimpse(df_CPUE_kg_day)

Yield for the period that matches Ecopath model, post 2000 annual average

df_yield_2000_2019 <- df_CPUE_kg_day %>% 
  filter(Year > 2000) %>% 
  group_by(Species) %>% 
  summarise(yield = sum(total_catch_kg)/19)

df_yield_2000_2019

Add yield in tonnes per km2, same as g m2 and it is consistent with biomass_observed

1.474341e+12 m^2 for model domain

1.474341e+12/1e+6 = 1474341 km^2

507511 kg of minke whale per year from 2000 - 2019

507.5/1474341 = 0.0003442216 tonnes Minke whale yield per km2 from 2000 - 2019

15619.24 kg of baleen whale per year from 2000 - 2019

15.6/1474341 = 1.0581e-05 tonnes of baleen whale per km2 from 2000 - 2019

obs_yield <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0003442216, 0, 0, 1.0581e-05) # to add a yield_observed column in species_params

biomass_cutoff <- c(0.1, 1, 1, 1, 40, 4500, 1, 1, 200000, 10000, 135000, 600000, 490000, 3650000, 2250000) # to add a biomass_cutoff column in species_params

Steady-state model with no fishing effort calibrated to observed biomass values from McCormack et al. 2020 that represent an average state of the food web from 2010-2020.

# params <-  readRDS("stage1_steady_vXX.rds") # Incorrect orca and leopard seal w_max values
so_params <- readRDS("params_16_March_2023.rds") # Updated w_max for orca and leopard seals

# params <- setParams(setFishing(params, initial_effort = 0.2))

# species_params(params)$yield_observed <- obs_yield
# species_params(params)$biomass_cutoff <- biomass_cutoff
# species_params(params)$biomass_cutoffLow <- biomass_cutoff
# species_params(params)$biomass_cutoffHigh <- species_params(params)$w_max

species_params(so_params)$yield_observed <- obs_yield
species_params(so_params)$biomass_cutoff <- biomass_cutoff
species_params(so_params)$biomass_cutoffLow <- biomass_cutoff
species_params(so_params)$biomass_cutoffHigh <- species_params(so_params)$w_max

# species_params(params) |> dplyr::select(species, yield_observed, w_mat, w_max, R_max)
species_params(so_params) |> dplyr::select(species, yield_observed, w_mat, w_max, R_max)

# groups |> dplyr::select(species, yield_observed, biomass_observed, biomass_cutoff, w_min)

# yield_observed <- c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.0003442216, 0, 0, 1.0581e-05) # to add a yield_observed column in species_params
# 
# species_params(params)$yield_observed <- yield_observed

# biomass_cutoff <- c(0.1, 1, 1, 1, 40, 4500, 1, 1, 200000, 10000, 135000, 600000, 490000, 3650000, 2250000) # to add a biomass_cutoff column in species_params
# 
# species_params(params)$biomass_cutoff <- biomass_cutoff

Update max and mat size for orca using sealifebase values rather than the previous w_max of 6t

# species_params(params)$species[13]
# species_params(params)$w_max[13]
# species_params(params)$w_mat[13]
# 
# species_params(params)$w_max[13] <- 10628034
# 
# species_params(params)$w_mat[13] <- 3198855
# 
# species_params(params)$w_max[13]
# species_params(params)$w_mat[13]
# species_params(params)$w_mat25[13]

Update max and mat size for leopard seal using sealifebase values. Max weight observed reported as 450kg. Using the sealifebase values for L-W conversion and max length results in an estimated max weight of 545875.2g, which is too high above the max weight observed.

# species_params(params)$species[9]
# species_params(params)$w_max[9]
# species_params(params)$w_mat[9]
# 
# species_params(params)$w_max[9] <- 450000
#  
# species_params(params)$w_max[9]
# species_params(params)$w_mat[9]
# species_params(params)$w_mat25[9]

Adjust w_mat values that were changed by default in newMultispeciesParams

# params_v1 <- params
# 
# params_v1@species_params$w_mat[params_v1@species_params$species == "large divers"]  <- params_v1@species_params$w_max[params_v1@species_params$species == "large divers"] * 0.9
# params_v1@species_params$w_mat[params_v1@species_params$species == "minke whales"]  <- params_v1@species_params$w_max[params_v1@species_params$species == "minke whales"] * 0.9
# params_v1@species_params$w_mat[params_v1@species_params$species == "orca"]  <- params_v1@species_params$w_max[params_v1@species_params$species == "orca"] * 0.9
# params_v1@species_params$w_mat[params_v1@species_params$species == "sperm whales"]  <- params_v1@species_params$w_max[params_v1@species_params$species == "sperm whales"] * 0.9
# 
# params_v1 <- setParams(params_v1)

# params_v2 <- params_v1
# 
# params_v2@species_params$w_min[params_v2@species_params$species == "small divers"]  <- params_v2@species_params$w_mat[params_v2@species_params$species == "small divers"] * 0.85
# params_v2@species_params$w_min[params_v2@species_params$species == "leopard seals"]  <- params_v2@species_params$w_mat[params_v2@species_params$species == "leopard seals"] * 0.85
# 
# params_v2 <- setParams(params_v2)

params <- steady(so_params)

Check gear params

gear_params(params)

Adjust catchability to only select fishing on species with catch data This is a only starting point

gear_params(params)$catchability <-  c(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0.49, 0, 0, 0.001)

Catchability: 1 value Effort changes through time (voyage frequency and length)

gear_params(params)$gear <-  c("Main", "Main", "Main", "Main", "Main", "Main",
                                   "Main", "Main", "Main", "Main", "Main", "Main",
                                   "Main", "Main", "Main")

gear_params(params)$l50 <-  c(1, 5, 5, 5, 10, 10, 10, 20, 20, 20, 20, 850, 600, 1500, 2200) # values for minke, orca, sperm and baleen are estimated based off `catch_lengths`, while all others are rough guesses, purely as the param won't work with values missing.

gear_params(params)$l25 <-  c(0.8, 4, 4, 4, 8, 8, 8, 16, 16, 16, 16, 800, 500, 1400, 1500)

gear_params(params)$sel_func <-  c("sigmoid_length", "sigmoid_length", "sigmoid_length", "sigmoid_length", "sigmoid_length", "sigmoid_length",
                                   "sigmoid_length", "sigmoid_length", "sigmoid_length", "sigmoid_length", "sigmoid_length", "sigmoid_length",
                                   "sigmoid_length", "sigmoid_length", "sigmoid_length")

gear_params(params)

params_v02 <- setParams(setFishing(params, initial_effort = 0.2))

params_v02 <- steady(params_v02)

sim_v01 <- project(params_v02, t_max = 100)
plot(sim_v01)
plotlyBiomass(sim_v01)
sim_v01@params@initial_effort

plotYieldVsSize(sim_v01, species = "baleen whales", catch = catch_lengths, 
                x_var = "Length")

# plotYieldVsSize(sim_v1, species = "sperm whales", catch = catch_lengths, 
#                 x_var = "Length")
# 
# plotYieldVsSize(sim_v1, species = "orca", catch = catch_lengths, 
#                 x_var = "Length")

plotYieldVsSize(sim_v01, species = "minke whales", catch = catch_lengths, 
                x_var = "Length")

getErrorCustom <- function(vary, params, dat, tol = 0.001, 
    timetorun = 10)
{
  params@species_params$R_max[1:15]<-10^vary[1:15] # R_max for 15 species
  params@species_params$erepro[1:15]<-vary[16:30] # erepro for 15 species
  params@species_params$interaction_resource[1:15] <- vary[31:45] # interaction_resource for 15 species
  params <- setParams(params)
  # interaction <- params@interaction
  # interaction[] <- matrix(vary[28:108],nrow = 9) # stop at 54 if looking only at 3 biggest species
  
  # params <- setInteraction(params,interaction)
    params <- projectToSteady(params, distance_func = distanceSSLogN, 
        tol = tol, t_max = 200, return_sim = F)
    
    sim <- project(params, t_max = timetorun, progress_bar = F)
    
    sim_biomass = rep(0, length(params@species_params$species))
    
        cutoffLow <- params@species_params$biomass_cutoffLow
    if (is.null(cutoffLow)) 
        cutoffLow <- rep(0, no_sp)
    cutoffLow[is.na(cutoffLow)] <- 0
    
        cutoffHigh <- params@species_params$biomass_cutoffHigh
    if (is.null(cutoffHigh)) 
        cutoffHigh <- rep(0, no_sp)
    cutoffHigh[is.na(cutoffHigh)] <- 0
        
    for (j in 1:length(sim_biomass)) {
        sim_biomass[j] = sum((sim@n[dim(sim@n)[1],j,] * params@w * 
            params@dw)[params@w >= cutoffLow[j] & cutoffHigh[j] >= params@w])
    }
    
     pred <- log(sim_biomass)
    dat <- log(dat)
    discrep <- pred - dat
    discrep <- (sum(discrep^2))
    return(discrep)
}
    

# create set of params for the optimisation process
tic()

params_optim <- params_v02

vary <- c(log10(params_optim@species_params$R_max),
          params_optim@species_params$erepro,
          params_optim@species_params$interaction_resource)

params_optim<-setParams(params_optim)

# set up workers
noCores <- parallel::detectCores() - 1 # keep some spare core
cl <- parallel::makeCluster(noCores, setup_timeout = 0.5)
setDefaultCluster(cl = cl)
clusterExport(cl, varlist = "cl",envir=environment())
clusterEvalQ(cl, {
  library(mizerExperimental)
  library(optimParallel)
})

optim_result <- optimParallel::optimParallel(par=vary,getErrorCustom,params=params_optim, 
                                             dat = params_optim@species_params$biomass_observed, method ="L-BFGS-B", 
                                             lower=c(rep(-15,15),rep(1e-7,15),rep(.1,15)),
                                             upper= c(rep(15,15),rep(1,15),rep(.99,15)),
                                             parallel=list(loginfo=TRUE, forward=TRUE))
stopCluster(cl)

toc()

saveRDS(optim_result, file="optim_result_v00.RDS")
# optim_result <- readRDS("optim_result_v00.RDS")

#put these new vals intospecies_params and go back to the top of this page to re-check the calibration 
species_params(params_optim)$R_max<-10^optim_result$par[1:15]
species_params(params_optim)$erepro<-optim_result$par[16:30]
species_params(params_optim)$interaction_resource <-optim_result$par[31:45]
 
sim_optim <- project(params_optim, t_max = 2000)
plotBiomass(sim_optim)

#put these new vals intospecies_params and go back to the top of this page to re-check the calibration 
species_params(params_optim)$R_max
species_params(params_optim)$erepro
species_params(params_optim)$interaction_resource

# this function adds a lower boundary to selected size
plotBiomassObservedVsModelCustom <- function (object, species = NULL, ratio = FALSE, log_scale = TRUE, 
    return_data = FALSE, labels = TRUE, show_unobserved = FALSE) 
{
    if (is(object, "MizerSim")) {
        params = object@params
        n <- finalN(object)
    }
    else if (is(object, "MizerParams")) {
        params = object
        n <- initialN(params)
    }
    else {
        stop("You have not provided a valid mizerSim or mizerParams object.")
    }
    sp_params <- params@species_params
    species = valid_species_arg(object, species)
    if (length(species) == 0) 
        stop("No species selected, please fix.")
    row_select = match(species, sp_params$species)
    if (!"biomass_observed" %in% names(sp_params)) {
        stop("You have not provided values for the column 'biomass_observed' ", 
            "in the mizerParams/mizerSim object.")
    }
    else if (!is.numeric(sp_params$biomass_observed)) {
        stop("The column 'biomass_observed' in the mizerParams/mizerSim object", 
            " is not numeric, please fix.")
    }
    else {
        biomass_observed = sp_params$biomass_observed
    }
    
    cutoffLow <- sp_params$biomass_cutoffLow[row_select]
    if (is.null(cutoffLow)) {
        cutoffLow = rep(0, length(species))
    }
    else if (!is.numeric(cutoffLow)) {
        stop("params@species_params$biomass_cutoffLow is not numeric, \",\n                 \"please fix.")
    }
    cutoffLow[is.na(cutoffLow)] <- 0
    
    cutoffHigh <- sp_params$biomass_cutoffHigh[row_select]
    if (is.null(cutoffHigh)) {
        cutoffHigh = rep(0, length(species))
    }
    else if (!is.numeric(cutoffHigh)) {
        stop("params@species_params$biomass_cutoffHigh is not numeric, \",\n                 \"please fix.")
    }
    cutoffHigh[is.na(cutoffHigh)] <- 0
    
    sim_biomass = rep(0, length(species))
    for (j in 1:length(species)) {
        sim_biomass[j] = sum((n[row_select[j], ] * params@w * 
            params@dw)[params@w >= cutoffLow[j] & cutoffHigh[j] >= params@w])
    }
    dummy = data.frame(species = species, model = sim_biomass, 
        observed = biomass_observed[row_select]) %>% mutate(species = factor(species, 
        levels = species), is_observed = !is.na(observed) & observed > 
        0, observed = case_when(is_observed ~ observed, !is_observed ~ 
        model), ratio = model/observed)
    if (sum(dummy$is_observed) == 0) {
        cat(paste("There are no observed biomasses to compare to model,", 
            "only plotting model biomasses.", sep = "\n"))
    }
    if (!show_unobserved) {
        dummy <- filter(dummy, is_observed)
    }
    if (return_data == TRUE) 
        return(dummy)
    tre <- round(sum(abs(1 - dummy$ratio)), digits = 3)
    caption <- paste0("Total relative error = ", tre)
    if (any(!dummy$is_observed)) {
        caption <- paste(caption, "\n Open circles represent species without biomass observation.")
    }
    if (ratio == FALSE) {
        gg <- ggplot(data = dummy, aes(x = observed, y = model, 
            colour = species, shape = is_observed)) + geom_abline(aes(intercept = 0, 
            slope = 1), colour = "purple", linetype = "dashed", 
            size = 1.3) + geom_point(size = 3) + labs(y = "model biomass [g]") + 
            coord_cartesian(ylim = range(dummy$model, dummy$observed))
    }
    else {
        gg <- ggplot(data = dummy, aes(x = observed, y = ratio, 
            colour = species, shape = is_observed)) + geom_hline(aes(yintercept = 1), 
            linetype = "dashed", colour = "purple", 
            size = 1.3) + geom_point(size = 3) + labs(y = "model biomass / observed biomass") + 
            coord_cartesian(ylim = range(dummy$ratio))
    }
    gg <- gg + labs(x = "observed biomass [g]", caption = caption) + 
        scale_colour_manual(values = getColours(params)[dummy$species]) + 
        scale_shape_manual(values = c(`TRUE` = 19, `FALSE` = 1)) + 
        guides(shape = "none")
    if (log_scale == TRUE & ratio == FALSE) {
        gg = gg + scale_x_log10() + scale_y_log10()
    }
    if (log_scale == TRUE & ratio == TRUE) {
        gg = gg + scale_x_log10()
    }
    if (labels == TRUE) {
        gg = gg + ggrepel::geom_label_repel(aes(label = species), 
            box.padding = 0.35, point.padding = 0.5, segment.color = "grey50", 
            show.legend = FALSE, max.overlaps = Inf, seed = 42)
    }
    gg
}

plotBiomassObservedVsModelCustom(sim_optim)

Use tuneParams() to investigate the species with erepro values too high

params_tuned_v01 <- tuneParams(params_optim)

params_tuned_v01@species_params$erepro

params_tuned_v01 <- steady(params_tuned_v01)

# create set of params for the optimisation process
tic()

params_optim <- params_tuned_v01

vary <- c(log10(params_optim@species_params$R_max),
          params_optim@species_params$erepro,
          params_optim@species_params$interaction_resource)

params_optim<-setParams(params_optim)

# set up workers
noCores <- parallel::detectCores() - 1 # keep some spare core
cl <- parallel::makeCluster(noCores, setup_timeout = 0.5)
setDefaultCluster(cl = cl)
clusterExport(cl, varlist = "cl",envir=environment())
clusterEvalQ(cl, {
  library(mizerExperimental)
  library(optimParallel)
})

optim_result <- optimParallel::optimParallel(par=vary,getErrorCustom,params=params_optim, 
                                             dat = params_optim@species_params$biomass_observed, method ="L-BFGS-B", 
                                             lower=c(rep(-15,15),rep(1e-7,15),rep(.1,15)),
                                             upper= c(rep(15,15),rep(0.99,15),rep(.99,15)),
                                             parallel=list(loginfo=TRUE, forward=TRUE))
stopCluster(cl)

toc()

saveRDS(optim_result, file="optim_result_v01.RDS")
# optim_result <- readRDS("optim_result_v00.RDS")

#put these new vals intospecies_params and go back to the top of this page to re-check the calibration 
species_params(params_optim)$R_max<-10^optim_result$par[1:15]
species_params(params_optim)$erepro<-optim_result$par[16:30]
species_params(params_optim)$interaction_resource <-optim_result$par[31:45]

species_params(params_optim)$R_max
species_params(params_optim)$erepro
species_params(params_optim)$interaction_resource

sim_optim_v02 <- project(params_optim, t_max = 1000)
plotBiomass(sim_optim_v02)

params_v03 <- steady(params_optim)

params_loop <- params_v03 |>
    matchBiomasses() |> steady() |> matchBiomasses() |> steady() |>
    matchBiomasses() |> steady() |> matchBiomasses() |> steady() |>
    matchBiomasses() |> steady() |> matchBiomasses() |> steady() |>
    matchBiomasses() |> steady() |> matchBiomasses() |> steady() |>
    matchBiomasses() |> steady() |> matchBiomasses() |> steady() |>
    matchBiomasses() |> steady() |> matchBiomasses() |> steady() |>
    matchBiomasses() |> steady() |> matchBiomasses() |> steady() |>
    matchBiomasses() |> steady() |> matchBiomasses() |> steady() |>
    matchBiomasses() |> steady() |> matchBiomasses() |> steady()

params_loop@species_params$erepro

sim_v02 <- project(params_loop, t_max = 2000)

plotBiomass(sim_v02)
plotBiomassObservedVsModelCustom(sim_v02)
plotBiomassObservedVsModel(sim_v02)
plotDiet(params_loop)

Save progress

# saveRDS(params_loop, "params_optim_v01.rds")

Try to ween the species of the resource spectra

params_loop@resource_params$w_pp_cutoff
params_loop@resource_params$kappa

resource_params(params_loop)[["w_pp_cutoff"]] <- 10

params_v04 <- setParams(params_loop)

params_v04@resource_params$w_pp_cutoff

sim_v03 <- project(params_v04, t_max = 200)

plotBiomass(sim_v03)
plotBiomassObservedVsModelCustom(sim_v03)
plotBiomassObservedVsModel(sim_v03)
plotDiet(params_v04)

getErrorCustom_v02 <- function(vary, params, dat, tol = 0.001, 
    timetorun = 10)
{
  params@species_params$R_max[1:15]<-10^vary[1:15] # R_max for 15 species
  params@species_params$erepro[1:15]<-vary[16:30] # erepro for 15 species
  params@species_params$interaction_resource[1:15] <- vary[31:45] # interaction_resource for 15 species
  params@resource_params$w_pp_cutoff[1] <- vary[46] # interaction_resource for 15 species
  params <- setParams(params)
  # interaction <- params@interaction
  # interaction[] <- matrix(vary[28:108],nrow = 9) # stop at 54 if looking only at 3 biggest species
  
  # params <- setInteraction(params,interaction)
    params <- projectToSteady(params, distance_func = distanceSSLogN, 
        tol = tol, t_max = 200, return_sim = F)
    
    sim <- project(params, t_max = timetorun, progress_bar = F)
    
    sim_biomass = rep(0, length(params@species_params$species))
    
        cutoffLow <- params@species_params$biomass_cutoffLow
    if (is.null(cutoffLow)) 
        cutoffLow <- rep(0, no_sp)
    cutoffLow[is.na(cutoffLow)] <- 0
    
        cutoffHigh <- params@species_params$biomass_cutoffHigh
    if (is.null(cutoffHigh)) 
        cutoffHigh <- rep(0, no_sp)
    cutoffHigh[is.na(cutoffHigh)] <- 0
        
    for (j in 1:length(sim_biomass)) {
        sim_biomass[j] = sum((sim@n[dim(sim@n)[1],j,] * params@w * 
            params@dw)[params@w >= cutoffLow[j] & cutoffHigh[j] >= params@w])
    }
    
     pred <- log(sim_biomass)
    dat <- log(dat)
    discrep <- pred - dat
    discrep <- (sum(discrep^2))
    return(discrep)
}
    

# create set of params for the optimisation process
tic()

params_optim <- params_v04

vary <- c(log10(params_optim@species_params$R_max),
          params_optim@species_params$erepro,
          params_optim@species_params$interaction_resource,
          params_optim@resource_params$w_pp_cutoff)

params_optim<-setParams(params_optim)

# set up workers
noCores <- parallel::detectCores() - 1 # keep some spare core
cl <- parallel::makeCluster(noCores, setup_timeout = 0.5)
setDefaultCluster(cl = cl)
clusterExport(cl, varlist = "cl",envir=environment())
clusterEvalQ(cl, {
  library(mizerExperimental)
  library(optimParallel)
})

optim_result <- optimParallel::optimParallel(par=vary,getErrorCustom_v02,params=params_optim, 
                                             dat = params_optim@species_params$biomass_observed, method ="L-BFGS-B", 
                                             lower=c(rep(-15,15),rep(1e-7,15),rep(.1,15), 0.1),
                                             upper= c(rep(15,15),rep(0.99,15),rep(.99,15),1000),
                                             parallel=list(loginfo=TRUE, forward=TRUE))
stopCluster(cl)

toc()

saveRDS(optim_result, file="optim_result_v02.RDS")
# optim_result <- readRDS("optim_result_v00.RDS")

#put these new vals intospecies_params and go back to the top of this page to re-check the calibration 
species_params(params_optim)$R_max<-10^optim_result$par[1:15]
species_params(params_optim)$erepro<-optim_result$par[16:30]
species_params(params_optim)$interaction_resource <-optim_result$par[31:45]
resource_params(params_optim)$w_pp_cutoff <- optim_result$par[46]

species_params(params_optim)$R_max
species_params(params_optim)$erepro
species_params(params_optim)$interaction_resource

sim_v04 <- project(params_optim, t_max = 1000)

plotBiomass(sim_v04)
plotBiomassObservedVsModelCustom(sim_v04)
plotBiomassObservedVsModel(sim_v04)
plotDiet(params_optim)

params_steady <- params_optim

initialN(params_steady) <- sim_v04@n[dim(sim_v04@n)[1],,]

plotDiet(params_steady)

sim_v05 <- project(params_steady, t_max = 100)

plotBiomass(sim_v05)

plotGrowthCurves(sim_v05, species_panel = T)

params_tuned_v02 <- tuneParams(params_steady)

params_tuned_v02 <- steady(params_tuned_v02)

params_tuned_v02@species_params$erepro

species w_min w_mat w_max beta k_vb h min_depth max_depth water.column.use p_time_prydz pc_annual_offspring 5 salps 3.162278e-05 0.2511886 25.11886 1000 NA 33 0 500 non DVM 1 NA biomass_observed 5 0.652

species_params <- data.frame(
    species = "salps",
    w_min = 3.162278e-05,
    w_max = 25.11886,
    w_mat = 0.2511886,
    beta = 1000000,
    sigma = 2,
    biomass_observed = 0.652,
    k_vb = 0.4,
    yield_observed = 0,
    biomass_cutoff = 0.1,
    biomass_cutoffLow = 0.1,
    biomass_cutoffHigh = 25.11886
)

params_v05 <- addSpecies(params_tuned_v02, species_params)
plotSpectra(params_v05)

params_v05@species_params$erepro

# function running tuneParams function in a row for a quick start to a calibration
fastCalib <- function(params, match = F)
{
params <- calibrateBiomass(params) # changes kappa and rmax
if(match) params <- matchBiomasses(params) # set rmax to inf and adjust erepro
params <- steady(params, tol = 0.001)
sim <- project(params, t_max = 1000)
return(sim)
}

sim_v06 <- fastCalib(params_v05)

plotBiomass(sim_v06)
plotDiet(params_v05)

so_theta <- as.matrix(params_v05@interaction)

salp_row <- c(1,
0.1250000,
0.1250000,
0.1373626,
0.0500000,
0.1000000,
0.2500000,
0.2500000,
0.1000000,
0.1666667,
0.1666667,
0.0500000,
0.2000000,
0.1250000,
0.2000000,
1)

so_theta[16,] # automated salps interaction when added using addSpecies()

so_theta[16,] <- salp_row

so_theta[,16] # automated salps interaction when added using addSpecies()

so_theta[,16] <- salp_row

so_theta

params_v06 <- setParams(params_v05, interaction = so_theta)

params_v06@interaction

sim_v07 <- fastCalib(params_v06)

plotBiomass(sim_v07)

initialN(params_v06) <- sim_v07@n[dim(sim_v07@n)[1],,]

sim_v08 <- project(params_v06, t_max = 100)
plotBiomass(sim_v08)
plotDiet(params_v06)

Fill in missing and adjust beta values for zooplankton groups using Heneghan et al. 2020 (10.1016/j.ecolmodel.2020.109265)

Need a value for microzooplankton microzooplankton (in McCormack et al. 2020) is composed of Heterotrophic dinoflagellates, tintinnids, ciliates, copepod nauplii Heneghan et al. log10PPMR values for: Hetero.Flagellates = 0.2–0.72 -> 0.46 Hetero.Ciliates = 2.5–2.9 -> 2.7 Mean: (2.7+0.46)/2 = 1.58 10^1.58 = 38.01894

Need to adjust values for mesozoo, other macrozoo, euphausiids, salps

Heneghan et al. log10PPMR values (midpoints for range) for: salps = 6.8–11.7 -> 9.25 10^9.25 = 1778279410

euphausiids = 6.6–7.8 -> 7.2 10^7.2 = 15848932

mesozoo (copepods) Omni.Cop. = 3.6–4.6 -> 4.1 Carn.Cop. = 0.8–1.9 -> 1.35 Mean: (4.1+1.35)/2 = 2.725 10^2.725 = 530.8844

other macrozoo () Chaetognaths = 1.9–3.4 -> 2.65 10^2.65 = 446.6836

params_tuned_v03 <- tuneParams(params_v06)
params_tuned_v04 <- tuneParams(params_tuned_v03) # updated salps PPMR to match empirical estimates from Heneghan et al 2020

getErrorCustom_v3 <- function(vary, params, dat, tol = 0.001, 
    timetorun = 10)
{
  params@species_params$R_max[1:16]<-10^vary[1:16] # R_max for 15 species
  params@species_params$erepro[1:16]<-vary[17:32] # erepro for 15 species
  params@species_params$interaction_resource[1:16] <- vary[33:48] # interaction_resource for 15 species
  params <- setParams(params)
  # interaction <- params@interaction
  # interaction[] <- matrix(vary[28:108],nrow = 9) # stop at 54 if looking only at 3 biggest species
  
  # params <- setInteraction(params,interaction)
    params <- projectToSteady(params, distance_func = distanceSSLogN, 
        tol = tol, t_max = 200, return_sim = F)
    
    sim <- project(params, t_max = timetorun, progress_bar = F)
    
    sim_biomass = rep(0, length(params@species_params$species))
    
        cutoffLow <- params@species_params$biomass_cutoffLow
    if (is.null(cutoffLow)) 
        cutoffLow <- rep(0, no_sp)
    cutoffLow[is.na(cutoffLow)] <- 0
    
        cutoffHigh <- params@species_params$biomass_cutoffHigh
    if (is.null(cutoffHigh)) 
        cutoffHigh <- rep(0, no_sp)
    cutoffHigh[is.na(cutoffHigh)] <- 0
        
    for (j in 1:length(sim_biomass)) {
        sim_biomass[j] = sum((sim@n[dim(sim@n)[1],j,] * params@w * 
            params@dw)[params@w >= cutoffLow[j] & cutoffHigh[j] >= params@w])
    }
    
     pred <- log(sim_biomass)
    dat <- log(dat)
    discrep <- pred - dat
    discrep <- (sum(discrep^2))
    return(discrep)
}

# create set of params for the optimisation process
tic()

params_optim <- params_tuned_v04

vary <- c(log10(params_optim@species_params$R_max),
          params_optim@species_params$erepro,
          params_optim@species_params$interaction_resource)

params_optim<-setParams(params_optim)

# set up workers
noCores <- parallel::detectCores() - 1 # keep some spare core
cl <- parallel::makeCluster(noCores, setup_timeout = 0.5)
setDefaultCluster(cl = cl)
clusterExport(cl, varlist = "cl",envir=environment())
clusterEvalQ(cl, {
  library(mizerExperimental)
  library(optimParallel)
})

optim_result <- optimParallel::optimParallel(par=vary,getErrorCustom_v3,params=params_optim, 
                                             dat = params_optim@species_params$biomass_observed, method ="L-BFGS-B", 
                                             lower=c(rep(-15,16),rep(1e-7,16),rep(.1,16)),
                                             upper= c(rep(15,16),rep(0.99,16),rep(.99,16)),
                                             parallel=list(loginfo=TRUE, forward=TRUE))
stopCluster(cl)

toc()

saveRDS(optim_result, file="optim_result_v03.RDS")
# optim_result <- readRDS("optim_result_v00.RDS")

#put these new vals intospecies_params and go back to the top of this page to re-check the calibration 
species_params(params_optim)$R_max<-10^optim_result$par[1:16]
species_params(params_optim)$erepro<-optim_result$par[17:32]
species_params(params_optim)$interaction_resource <-optim_result$par[33:48]

species_params(params_optim)$R_max
species_params(params_optim)$erepro
species_params(params_optim)$interaction_resource

sim_v09 <- project(params_optim, t_max = 200)
plotBiomass(sim_v09)

params_v07 <- params_optim
initialN(params_v07) <- sim_v09@n[dim(sim_v09@n)[1],,]

sim_v10 <- project(params_v07, t_max = 200)
plotBiomass(sim_v10)
plotDiet(params_v07)

initialN(params_v07) <- sim_v10@n[dim(sim_v10@n)[1],,]

sim_v11 <- project(params_v07, t_max = 200)
plotBiomass(sim_v11)

params_tuned_v05 <- tuneParams(params_v07)

params_tuned_v05 <- steady(params_tuned_v05)

params_tuned_v05@species_params$erepro

# saveRDS(params_tuned_v05, "params_optim_v02.rds")

params_tuned_v05 <- readRDS("params_optim_v02.rds")

sim_v12 <- project(params_tuned_v05, t_max = 500)

initialN(params_tuned_v05) <- sim_v12@n[dim(sim_v12@n)[1],,]

sim_v13 <- project(params_tuned_v05, t_max = 500)

plotlyBiomass(sim_v13)
plotBiomassObservedVsModel(sim_v13)
plotBiomassObservedVsModelCustom(sim_v13)

Gradually reducing resource maximum size in tuneParams(), as doing it directly in the setParams() route will tell you the w_pp_cutoff has changed, but it will still incorporate a background resource up to the w_pp_cutoff that was originally used in newMultispeciesParams()

params_tuned_v06 <- tuneParams(params_tuned_v05)

params_tuned_v06 <- steady(params_tuned_v06)

params_tuned_v06@species_params$erepro
params_tuned_v06@species_params$R_max

params_tuned_v10 <- tuneParams(params_tuned_v06)

params_tuned_v10 <- steady(params_tuned_v10)

params_tuned_v10@species_params$erepro
params_tuned_v10@species_params$R_max

params_tuned_v10@resource_params$w_pp_cutoff

# saveRDS(params_tuned_v06, "params_optim_v02_w_pp_100.rds")
# saveRDS(params_tuned_v10, "params_optim_v03_w_pp_100.rds")

# create set of params for the optimisation process
tic()

params_optim <- params_tuned_v10

vary <- c(log10(params_optim@species_params$R_max),
          params_optim@species_params$erepro,
          params_optim@species_params$interaction_resource)

params_optim<-setParams(params_optim)

# set up workers
noCores <- parallel::detectCores() - 1 # keep some spare core
cl <- parallel::makeCluster(noCores, setup_timeout = 0.5)
setDefaultCluster(cl = cl)
clusterExport(cl, varlist = "cl",envir=environment())
clusterEvalQ(cl, {
  library(mizerExperimental)
  library(optimParallel)
})

optim_result <- optimParallel::optimParallel(par=vary,getErrorCustom_v3,params=params_optim,
                                             dat = params_optim@species_params$biomass_observed, method ="L-BFGS-B",
                                             lower=c(rep(-15,16),rep(1e-7,16),rep(.1,16)),
                                             upper= c(rep(15,16),rep(0.99,16),rep(.99,16)),
                                             parallel=list(loginfo=TRUE, forward=TRUE))
stopCluster(cl)

toc()

saveRDS(optim_result, file="optim_result_v04.RDS")

# optim_result <- readRDS("optim_result_v04.RDS")

#put these new vals intospecies_params and go back to the top of this page to re-check the calibration 
species_params(params_optim)$R_max<-10^optim_result$par[1:16]
species_params(params_optim)$erepro<-optim_result$par[17:32]
species_params(params_optim)$interaction_resource <-optim_result$par[33:48]

species_params(params_optim)$R_max
species_params(params_optim)$erepro
species_params(params_optim)$interaction_resource

sim_v14 <- project(params_optim, t_max = 2000)
plotBiomass(sim_v14)

params_v08 <- params_optim
initialN(params_v08) <- sim_v14@n[dim(sim_v14@n)[1],,]

# saveRDS(params_tuned_v07, file="params_optim_v04_w_pp_100.rds")
params_tuned_v07 <- readRDS("params_optim_v04_w_pp_100.rds")

Warning: cannot open compressed file 'params_optim_v04_w_pp_100.rds', probable reason 'No such file or directory'Error in gzfile(file, "rb") : cannot open the connection

sim_v15 <- project(params_tuned_v07, t_max = 2000)

[--------------------------------------------]   0% ETA:  1m
[--------------------------------------------]   0% ETA:  2m
[--------------------------------------------]   0% ETA:  1m
[--------------------------------------------]   0% ETA:  2m
[--------------------------------------------]   1% ETA:  2m
[--------------------------------------------]   1% ETA:  1m
[--------------------------------------------]   1% ETA:  2m
[>-------------------------------------------]   1% ETA:  2m
[>-------------------------------------------]   1% ETA:  1m
[>-------------------------------------------]   1% ETA:  2m
[>-------------------------------------------]   2% ETA:  2m
[>-------------------------------------------]   2% ETA:  1m
[>-------------------------------------------]   2% ETA:  2m
[>-------------------------------------------]   2% ETA:  1m
[>-------------------------------------------]   2% ETA:  2m
[>-------------------------------------------]   2% ETA:  1m
[>-------------------------------------------]   2% ETA:  2m
[>-------------------------------------------]   2% ETA:  1m
[>-------------------------------------------]   3% ETA:  2m
[>-------------------------------------------]   3% ETA:  1m
[>-------------------------------------------]   3% ETA:  2m
[>-------------------------------------------]   3% ETA:  1m
[>-------------------------------------------]   3% ETA:  2m
[>-------------------------------------------]   3% ETA:  1m
[>-------------------------------------------]   3% ETA:  2m
[>-------------------------------------------]   3% ETA:  1m
[>-------------------------------------------]   3% ETA:  2m
[>-------------------------------------------]   3% ETA:  1m
[=>------------------------------------------]   3% ETA:  1m
[=>------------------------------------------]   3% ETA:  2m
[=>------------------------------------------]   4% ETA:  1m
[=>------------------------------------------]   5% ETA:  1m
[=>------------------------------------------]   6% ETA:  1m
[==>-----------------------------------------]   6% ETA:  1m
[==>-----------------------------------------]   7% ETA:  1m
[==>-----------------------------------------]   8% ETA:  1m
[===>----------------------------------------]   8% ETA:  1m
[===>----------------------------------------]   9% ETA:  1m
[===>----------------------------------------]  10% ETA:  1m
[====>---------------------------------------]  10% ETA:  1m
[====>---------------------------------------]  11% ETA:  1m
[====>---------------------------------------]  12% ETA:  1m
[=====>--------------------------------------]  13% ETA:  1m
[=====>--------------------------------------]  14% ETA:  1m
[=====>--------------------------------------]  15% ETA:  1m
[======>-------------------------------------]  15% ETA:  1m
[======>-------------------------------------]  16% ETA:  1m
[======>-------------------------------------]  17% ETA:  1m
[=======>------------------------------------]  17% ETA:  1m
[=======>------------------------------------]  18% ETA:  1m
[=======>------------------------------------]  19% ETA:  1m
[========>-----------------------------------]  19% ETA:  1m
[========>-----------------------------------]  20% ETA:  1m
[========>-----------------------------------]  21% ETA:  1m
[========>-----------------------------------]  22% ETA:  1m
[=========>----------------------------------]  22% ETA:  1m
[=========>----------------------------------]  23% ETA:  1m
[=========>----------------------------------]  24% ETA:  1m
[==========>---------------------------------]  24% ETA:  1m
[==========>---------------------------------]  25% ETA:  1m
[==========>---------------------------------]  26% ETA:  1m
[===========>--------------------------------]  26% ETA:  1m
[===========>--------------------------------]  27% ETA:  1m
[===========>--------------------------------]  28% ETA:  1m
[============>-------------------------------]  28% ETA:  1m
[============>-------------------------------]  29% ETA:  1m
[============>-------------------------------]  30% ETA:  1m
[============>-------------------------------]  31% ETA:  1m
[=============>------------------------------]  31% ETA:  1m
[=============>------------------------------]  32% ETA:  1m
[=============>------------------------------]  33% ETA:  1m
[==============>-----------------------------]  33% ETA:  1m
[==============>-----------------------------]  34% ETA:  1m
[==============>-----------------------------]  35% ETA:  1m
[===============>----------------------------]  35% ETA:  1m
[===============>----------------------------]  36% ETA:  1m
[===============>----------------------------]  37% ETA:  1m
[================>---------------------------]  38% ETA:  1m
[================>---------------------------]  39% ETA:  1m
[================>---------------------------]  40% ETA:  1m
[=================>--------------------------]  40% ETA:  1m
[=================>--------------------------]  41% ETA:  1m
[=================>--------------------------]  42% ETA:  1m
[==================>-------------------------]  42% ETA:  1m
[==================>-------------------------]  42% ETA: 50s
[==================>-------------------------]  43% ETA: 50s
[==================>-------------------------]  43% ETA: 49s
[==================>-------------------------]  43% ETA: 50s
[==================>-------------------------]  43% ETA: 49s
[==================>-------------------------]  44% ETA: 49s
[==================>-------------------------]  44% ETA: 48s
[==================>-------------------------]  44% ETA: 49s
[==================>-------------------------]  44% ETA: 48s
[===================>------------------------]  44% ETA: 48s
[===================>------------------------]  45% ETA: 48s
[===================>------------------------]  45% ETA: 47s
[===================>------------------------]  46% ETA: 47s
[===================>------------------------]  46% ETA: 46s
[===================>------------------------]  47% ETA: 46s
[====================>-----------------------]  47% ETA: 46s
[====================>-----------------------]  47% ETA: 45s
[====================>-----------------------]  48% ETA: 45s
[====================>-----------------------]  48% ETA: 44s
[====================>-----------------------]  49% ETA: 44s
[=====================>----------------------]  49% ETA: 44s
[=====================>----------------------]  49% ETA: 43s
[=====================>----------------------]  50% ETA: 44s
[=====================>----------------------]  50% ETA: 43s
[=====================>----------------------]  51% ETA: 43s
[=====================>----------------------]  51% ETA: 42s
[======================>---------------------]  51% ETA: 42s
[======================>---------------------]  52% ETA: 42s
[======================>---------------------]  52% ETA: 41s
[======================>---------------------]  53% ETA: 41s
[======================>---------------------]  53% ETA: 40s
[======================>---------------------]  53% ETA: 41s
[======================>---------------------]  53% ETA: 40s
[=======================>--------------------]  53% ETA: 40s
[=======================>--------------------]  54% ETA: 40s
[=======================>--------------------]  54% ETA: 39s
[=======================>--------------------]  55% ETA: 39s
[=======================>--------------------]  55% ETA: 38s
[=======================>--------------------]  56% ETA: 38s
[========================>-------------------]  56% ETA: 38s
[========================>-------------------]  56% ETA: 37s
[========================>-------------------]  56% ETA: 38s
[========================>-------------------]  56% ETA: 37s
[========================>-------------------]  57% ETA: 37s
[========================>-------------------]  57% ETA: 36s
[========================>-------------------]  58% ETA: 36s
[=========================>------------------]  58% ETA: 36s
[=========================>------------------]  59% ETA: 36s
[=========================>------------------]  59% ETA: 35s
[=========================>------------------]  60% ETA: 35s
[=========================>------------------]  60% ETA: 34s
[==========================>-----------------]  60% ETA: 34s
[==========================>-----------------]  61% ETA: 34s
[==========================>-----------------]  61% ETA: 33s
[==========================>-----------------]  62% ETA: 33s
[==========================>-----------------]  62% ETA: 32s
[===========================>----------------]  63% ETA: 32s
[===========================>----------------]  63% ETA: 31s
[===========================>----------------]  64% ETA: 31s
[===========================>----------------]  64% ETA: 30s
[===========================>----------------]  65% ETA: 30s
[============================>---------------]  65% ETA: 30s
[============================>---------------]  65% ETA: 29s
[============================>---------------]  66% ETA: 29s
[============================>---------------]  67% ETA: 29s
[============================>---------------]  67% ETA: 28s
[=============================>--------------]  67% ETA: 28s
[=============================>--------------]  68% ETA: 28s
[=============================>--------------]  68% ETA: 27s
[=============================>--------------]  69% ETA: 27s
[=============================>--------------]  69% ETA: 26s
[==============================>-------------]  69% ETA: 26s
[==============================>-------------]  70% ETA: 26s
[==============================>-------------]  70% ETA: 25s
[==============================>-------------]  71% ETA: 25s
[==============================>-------------]  71% ETA: 24s
[==============================>-------------]  72% ETA: 24s
[===============================>------------]  72% ETA: 24s
[===============================>------------]  72% ETA: 23s
[===============================>------------]  73% ETA: 23s
[===============================>------------]  74% ETA: 23s
[===============================>------------]  74% ETA: 22s
[================================>-----------]  74% ETA: 22s
[================================>-----------]  75% ETA: 22s
[================================>-----------]  75% ETA: 21s
[================================>-----------]  76% ETA: 21s
[================================>-----------]  76% ETA: 20s
[=================================>----------]  76% ETA: 20s
[=================================>----------]  77% ETA: 20s
[=================================>----------]  77% ETA: 19s
[=================================>----------]  78% ETA: 19s
[=================================>----------]  78% ETA: 18s
[==================================>---------]  78% ETA: 18s
[==================================>---------]  79% ETA: 18s
[==================================>---------]  79% ETA: 17s
[==================================>---------]  80% ETA: 17s
[==================================>---------]  81% ETA: 17s
[==================================>---------]  81% ETA: 16s
[===================================>--------]  81% ETA: 16s
[===================================>--------]  82% ETA: 16s
[===================================>--------]  82% ETA: 15s
[===================================>--------]  83% ETA: 15s
[====================================>-------]  83% ETA: 15s
[====================================>-------]  83% ETA: 14s
[====================================>-------]  84% ETA: 14s
[====================================>-------]  84% ETA: 13s
[====================================>-------]  85% ETA: 13s
[=====================================>------]  85% ETA: 13s
[=====================================>------]  85% ETA: 12s
[=====================================>------]  86% ETA: 12s
[=====================================>------]  87% ETA: 12s
[=====================================>------]  87% ETA: 11s
[======================================>-----]  88% ETA: 11s
[======================================>-----]  88% ETA: 10s
[======================================>-----]  89% ETA: 10s
[======================================>-----]  89% ETA:  9s
[======================================>-----]  90% ETA:  9s
[=======================================>----]  90% ETA:  9s
[=======================================>----]  90% ETA:  8s
[=======================================>----]  91% ETA:  8s
[=======================================>----]  91% ETA:  7s
[=======================================>----]  92% ETA:  7s
[========================================>---]  92% ETA:  7s
[========================================>---]  92% ETA:  6s
[========================================>---]  93% ETA:  6s
[========================================>---]  94% ETA:  6s
[========================================>---]  94% ETA:  5s
[=========================================>--]  94% ETA:  5s
[=========================================>--]  95% ETA:  5s
[=========================================>--]  95% ETA:  4s
[=========================================>--]  96% ETA:  4s
[=========================================>--]  96% ETA:  3s
[=========================================>--]  97% ETA:  3s
[==========================================>-]  97% ETA:  3s
[==========================================>-]  97% ETA:  2s
[==========================================>-]  98% ETA:  2s
[==========================================>-]  98% ETA:  1s
[==========================================>-]  99% ETA:  1s
[===========================================>]  99% ETA:  1s
[===========================================>]  99% ETA:  0s
[===========================================>] 100% ETA:  0s

plotBiomass(sim_v15)

plotDiet(params_tuned_v07)

initialN(params_tuned_v07) <- sim_v15@n[dim(sim_v15@n)[1],,]

sim_v16 <- project(params_tuned_v07, t_max = 500)

[>-------------------------------------------]   3% ETA:  7s
[=>------------------------------------------]   4% ETA:  7s
[=>------------------------------------------]   4% ETA: 11s
[=>------------------------------------------]   5% ETA: 10s
[=>------------------------------------------]   6% ETA: 10s
[==>-----------------------------------------]   6% ETA: 10s
[==>-----------------------------------------]   6% ETA:  9s
[==>-----------------------------------------]   6% ETA: 10s
[==>-----------------------------------------]   7% ETA:  9s
[==>-----------------------------------------]   8% ETA:  9s
[===>----------------------------------------]   8% ETA:  9s
[===>----------------------------------------]   9% ETA:  9s
[===>----------------------------------------]   9% ETA:  8s
[===>----------------------------------------]  10% ETA:  8s
[====>---------------------------------------]  10% ETA:  8s
[====>---------------------------------------]  11% ETA:  8s
[====>---------------------------------------]  12% ETA:  8s
[=====>--------------------------------------]  13% ETA:  8s
[=====>--------------------------------------]  13% ETA:  7s
[=====>--------------------------------------]  14% ETA:  7s
[=====>--------------------------------------]  15% ETA:  7s
[======>-------------------------------------]  15% ETA:  7s
[======>-------------------------------------]  16% ETA:  7s
[======>-------------------------------------]  17% ETA:  7s
[=======>------------------------------------]  17% ETA:  7s
[=======>------------------------------------]  18% ETA:  7s
[=======>------------------------------------]  19% ETA:  7s
[========>-----------------------------------]  19% ETA:  7s
[========>-----------------------------------]  20% ETA:  6s
[========>-----------------------------------]  21% ETA:  6s
[========>-----------------------------------]  22% ETA:  6s
[=========>----------------------------------]  22% ETA:  6s
[=========>----------------------------------]  23% ETA:  6s
[=========>----------------------------------]  24% ETA:  6s
[==========>---------------------------------]  24% ETA:  6s
[==========>---------------------------------]  25% ETA:  6s
[==========>---------------------------------]  26% ETA:  6s
[===========>--------------------------------]  26% ETA:  6s
[===========>--------------------------------]  27% ETA:  6s
[===========>--------------------------------]  28% ETA:  6s
[============>-------------------------------]  29% ETA:  5s
[============>-------------------------------]  30% ETA:  5s
[============>-------------------------------]  31% ETA:  5s
[=============>------------------------------]  31% ETA:  5s
[=============>------------------------------]  32% ETA:  5s
[=============>------------------------------]  33% ETA:  5s
[==============>-----------------------------]  33% ETA:  5s
[==============>-----------------------------]  34% ETA:  5s
[==============>-----------------------------]  35% ETA:  5s
[===============>----------------------------]  35% ETA:  5s
[===============>----------------------------]  36% ETA:  5s
[===============>----------------------------]  37% ETA:  5s
[================>---------------------------]  38% ETA:  5s
[================>---------------------------]  39% ETA:  4s
[================>---------------------------]  40% ETA:  4s
[=================>--------------------------]  40% ETA:  4s
[=================>--------------------------]  41% ETA:  4s
[=================>--------------------------]  42% ETA:  4s
[==================>-------------------------]  42% ETA:  4s
[==================>-------------------------]  43% ETA:  4s
[==================>-------------------------]  44% ETA:  4s
[===================>------------------------]  45% ETA:  4s
[===================>------------------------]  46% ETA:  4s
[===================>------------------------]  47% ETA:  4s
[====================>-----------------------]  47% ETA:  4s
[====================>-----------------------]  48% ETA:  4s
[====================>-----------------------]  49% ETA:  4s
[=====================>----------------------]  49% ETA:  4s
[=====================>----------------------]  50% ETA:  4s
[=====================>----------------------]  51% ETA:  4s
[======================>---------------------]  51% ETA:  4s
[======================>---------------------]  52% ETA:  4s
[======================>---------------------]  52% ETA:  3s
[======================>---------------------]  53% ETA:  3s
[=======================>--------------------]  53% ETA:  3s
[=======================>--------------------]  54% ETA:  3s
[=======================>--------------------]  55% ETA:  3s
[========================>-------------------]  56% ETA:  3s
[========================>-------------------]  57% ETA:  3s
[========================>-------------------]  58% ETA:  3s
[=========================>------------------]  58% ETA:  3s
[=========================>------------------]  59% ETA:  3s
[=========================>------------------]  60% ETA:  3s
[==========================>-----------------]  60% ETA:  3s
[==========================>-----------------]  61% ETA:  3s
[==========================>-----------------]  62% ETA:  3s
[===========================>----------------]  63% ETA:  3s
[===========================>----------------]  64% ETA:  3s
[===========================>----------------]  65% ETA:  2s
[============================>---------------]  65% ETA:  2s
[============================>---------------]  66% ETA:  2s
[============================>---------------]  67% ETA:  2s
[=============================>--------------]  67% ETA:  2s
[=============================>--------------]  68% ETA:  2s
[=============================>--------------]  69% ETA:  2s
[==============================>-------------]  69% ETA:  2s
[==============================>-------------]  70% ETA:  2s
[==============================>-------------]  71% ETA:  2s
[===============================>------------]  72% ETA:  2s
[===============================>------------]  73% ETA:  2s
[===============================>------------]  74% ETA:  2s
[================================>-----------]  74% ETA:  2s
[================================>-----------]  75% ETA:  2s
[================================>-----------]  76% ETA:  2s
[=================================>----------]  76% ETA:  2s
[=================================>----------]  77% ETA:  2s
[=================================>----------]  78% ETA:  2s
[=================================>----------]  78% ETA:  1s
[==================================>---------]  78% ETA:  1s
[==================================>---------]  79% ETA:  1s
[==================================>---------]  80% ETA:  1s
[==================================>---------]  81% ETA:  1s
[===================================>--------]  81% ETA:  1s
[===================================>--------]  82% ETA:  1s
[===================================>--------]  83% ETA:  1s
[====================================>-------]  83% ETA:  1s
[====================================>-------]  84% ETA:  1s
[====================================>-------]  85% ETA:  1s
[=====================================>------]  85% ETA:  1s
[=====================================>------]  86% ETA:  1s
[=====================================>------]  87% ETA:  1s
[======================================>-----]  88% ETA:  1s
[======================================>-----]  89% ETA:  1s
[======================================>-----]  90% ETA:  1s
[=======================================>----]  90% ETA:  1s
[=======================================>----]  91% ETA:  1s
[=======================================>----]  92% ETA:  1s
[========================================>---]  92% ETA:  1s
[========================================>---]  93% ETA:  1s
[========================================>---]  93% ETA:  0s
[========================================>---]  94% ETA:  0s
[=========================================>--]  94% ETA:  0s
[=========================================>--]  95% ETA:  0s
[=========================================>--]  96% ETA:  0s
[==========================================>-]  97% ETA:  0s
[==========================================>-]  98% ETA:  0s
[==========================================>-]  99% ETA:  0s
[===========================================>]  99% ETA:  0s
[===========================================>] 100% ETA:  0s

plotlyBiomass(sim_v16)

plotBiomassObservedVsModel(sim_v16)

plotlyBiomass(params_loop)

Error: is(object = sim, class2 = "MizerSim") is not TRUE

sim_v17 <- project(params_loop, t_max = 500)

[>-------------------------------------------]   3% ETA:  6s
[=>------------------------------------------]   4% ETA:  6s
[=>------------------------------------------]   4% ETA:  9s
[=>------------------------------------------]   5% ETA:  9s
[=>------------------------------------------]   6% ETA:  9s
[==>-----------------------------------------]   6% ETA:  8s
[==>-----------------------------------------]   7% ETA:  8s
[==>-----------------------------------------]   7% ETA: 10s
[==>-----------------------------------------]   7% ETA:  9s
[==>-----------------------------------------]   8% ETA:  9s
[===>----------------------------------------]   8% ETA:  9s
[===>----------------------------------------]   9% ETA:  9s
[===>----------------------------------------]  10% ETA:  9s
[===>----------------------------------------]  10% ETA:  8s
[====>---------------------------------------]  10% ETA:  8s
[====>---------------------------------------]  11% ETA:  8s
[====>---------------------------------------]  12% ETA:  8s
[=====>--------------------------------------]  13% ETA:  8s
[=====>--------------------------------------]  14% ETA:  8s
[=====>--------------------------------------]  14% ETA:  7s
[=====>--------------------------------------]  15% ETA:  7s
[======>-------------------------------------]  15% ETA:  7s
[======>-------------------------------------]  16% ETA:  7s
[======>-------------------------------------]  17% ETA:  7s
[=======>------------------------------------]  17% ETA:  7s
[=======>------------------------------------]  18% ETA:  7s
[=======>------------------------------------]  19% ETA:  7s
[========>-----------------------------------]  19% ETA:  7s
[========>-----------------------------------]  20% ETA:  7s
[========>-----------------------------------]  20% ETA:  6s
[========>-----------------------------------]  21% ETA:  6s
[========>-----------------------------------]  22% ETA:  6s
[=========>----------------------------------]  22% ETA:  6s
[=========>----------------------------------]  23% ETA:  6s
[=========>----------------------------------]  24% ETA:  6s
[==========>---------------------------------]  24% ETA:  6s
[==========>---------------------------------]  25% ETA:  6s
[==========>---------------------------------]  26% ETA:  6s
[===========>--------------------------------]  26% ETA:  6s
[===========>--------------------------------]  27% ETA:  6s
[===========>--------------------------------]  28% ETA:  6s
[============>-------------------------------]  29% ETA:  6s
[============>-------------------------------]  29% ETA:  5s
[============>-------------------------------]  30% ETA:  5s
[============>-------------------------------]  31% ETA:  5s
[=============>------------------------------]  31% ETA:  5s
[=============>------------------------------]  32% ETA:  5s
[=============>------------------------------]  33% ETA:  5s
[==============>-----------------------------]  33% ETA:  5s
[==============>-----------------------------]  34% ETA:  5s
[==============>-----------------------------]  35% ETA:  5s
[===============>----------------------------]  35% ETA:  5s
[===============>----------------------------]  36% ETA:  5s
[===============>----------------------------]  37% ETA:  5s
[================>---------------------------]  38% ETA:  5s
[================>---------------------------]  39% ETA:  5s
[================>---------------------------]  40% ETA:  5s
[=================>--------------------------]  40% ETA:  5s
[=================>--------------------------]  41% ETA:  5s
[=================>--------------------------]  42% ETA:  5s
[==================>-------------------------]  42% ETA:  5s
[==================>-------------------------]  43% ETA:  5s
[==================>-------------------------]  43% ETA:  4s
[==================>-------------------------]  44% ETA:  4s
[===================>------------------------]  45% ETA:  4s
[===================>------------------------]  46% ETA:  4s
[===================>------------------------]  47% ETA:  4s
[====================>-----------------------]  47% ETA:  4s
[====================>-----------------------]  48% ETA:  4s
[====================>-----------------------]  49% ETA:  4s
[=====================>----------------------]  49% ETA:  4s
[=====================>----------------------]  50% ETA:  4s
[=====================>----------------------]  51% ETA:  4s
[======================>---------------------]  51% ETA:  4s
[======================>---------------------]  52% ETA:  4s
[======================>---------------------]  53% ETA:  4s
[=======================>--------------------]  53% ETA:  4s
[=======================>--------------------]  54% ETA:  4s
[=======================>--------------------]  55% ETA:  4s
[========================>-------------------]  56% ETA:  4s
[========================>-------------------]  57% ETA:  4s
[========================>-------------------]  57% ETA:  3s
[========================>-------------------]  58% ETA:  3s
[=========================>------------------]  58% ETA:  3s
[=========================>------------------]  59% ETA:  3s
[=========================>------------------]  60% ETA:  3s
[==========================>-----------------]  60% ETA:  3s
[==========================>-----------------]  61% ETA:  3s
[==========================>-----------------]  62% ETA:  3s
[===========================>----------------]  63% ETA:  3s
[===========================>----------------]  64% ETA:  3s
[===========================>----------------]  65% ETA:  3s
[============================>---------------]  65% ETA:  3s
[============================>---------------]  66% ETA:  3s
[============================>---------------]  67% ETA:  3s
[=============================>--------------]  67% ETA:  3s
[=============================>--------------]  68% ETA:  3s
[=============================>--------------]  69% ETA:  3s
[==============================>-------------]  69% ETA:  3s
[==============================>-------------]  70% ETA:  2s
[==============================>-------------]  71% ETA:  2s
[===============================>------------]  72% ETA:  2s
[===============================>------------]  73% ETA:  2s
[===============================>------------]  74% ETA:  2s
[================================>-----------]  74% ETA:  2s
[================================>-----------]  75% ETA:  2s
[================================>-----------]  76% ETA:  2s
[=================================>----------]  76% ETA:  2s
[=================================>----------]  77% ETA:  2s
[=================================>----------]  78% ETA:  2s
[==================================>---------]  78% ETA:  2s
[==================================>---------]  79% ETA:  2s
[==================================>---------]  80% ETA:  2s
[==================================>---------]  81% ETA:  2s
[===================================>--------]  81% ETA:  2s
[===================================>--------]  81% ETA:  1s
[===================================>--------]  82% ETA:  1s
[===================================>--------]  83% ETA:  1s
[====================================>-------]  83% ETA:  1s
[====================================>-------]  84% ETA:  1s
[====================================>-------]  85% ETA:  1s
[=====================================>------]  85% ETA:  1s
[=====================================>------]  86% ETA:  1s
[=====================================>------]  87% ETA:  1s
[======================================>-----]  88% ETA:  1s
[======================================>-----]  89% ETA:  1s
[======================================>-----]  90% ETA:  1s
[=======================================>----]  90% ETA:  1s
[=======================================>----]  91% ETA:  1s
[=======================================>----]  92% ETA:  1s
[========================================>---]  92% ETA:  1s
[========================================>---]  93% ETA:  1s
[========================================>---]  94% ETA:  1s
[========================================>---]  94% ETA:  0s
[=========================================>--]  94% ETA:  0s
[=========================================>--]  95% ETA:  0s
[=========================================>--]  96% ETA:  0s
[==========================================>-]  97% ETA:  0s
[==========================================>-]  98% ETA:  0s
[==========================================>-]  99% ETA:  0s
[===========================================>]  99% ETA:  0s
[===========================================>] 100% ETA:  0s

plotlyBiomass(sim_v17)

plotBiomassObservedVsModel(sim_v17)

plotBiomassObservedVsModelCustom(sim_v17)

Error in plotBiomassObservedVsModelCustom(sim_v17) : 
  could not find function "plotBiomassObservedVsModelCustom"

params_loop <- steady(params_loop)

Convergence was achieved in 1.5 years.

Reduce interaction with resource gradually to see if that can increase predation on euphausiids

box.params <- params_loop

box.params@species_params$ppmr_min[box.params@species_params$species == "baleen whales"]  <- 4e6
box.params@species_params$ppmr_max[box.params@species_params$species == "baleen whales"] <-5e6
# box.params@species_params$pred_kernel_type[box.params@species_params$species == "baleen whales"] <- "box"

params_v15 <- tuneParams(box.params)

Listening on http://127.0.0.1:7954

NA